%% Simplex Dynamics
% Thomas Stevens, Lim Lab, opened 11/26/12

%Functions needed:
%quiver/quiver3
%gradient
%ndgrid
%TriScatteredInterp

%% First pass test

[X,Y] = meshgrid(1:1:5,1:1:5);
N = X+Y;
Xf = X./N;
Yf = Y./N;
%Payoff matrix
a = 2;
b = 1;
c = 2;
d = 1;

dX = a.*X + b.*Y;
dY = c.*X + d.*Y;
dN = dX+dY;
dXf = (N.*dX - X.*dN)./(N.^2);
dYf = (N.*dY - Y.*dN)./(N.^2);
%dXf + dYf should be 0.

%quiver(Xf,Yf,dXf,dYf);

%% 3-state model
grid_natural = [0 logspace(0,2,10)];
gridlen = length(grid_natural);
%griddim = linspace(0,9,10);
[X,Y,Z] = meshgrid(grid_natural,grid_natural,grid_natural);
N = X+Y+Z;
Xf = X./N;
Yf = Y./N;
Zf = Z./N;
%Paramters (Payoff matrix)
ka = 1;
%kd = 1;
fxx = 7;
fyy = 6;
fzz = 3;
    %x      %y    %z
A = [fxx     -ka     -ka;   %x
     -ka       fyy   -ka;   %y
     ka      ka     fzz];   %z

%% Linear payoffs = unique vector field in R2
dXYZ = A*[X(:)';Y(:)';Z(:)'];
dX = reshape(dXYZ(1,:),gridlen,gridlen,gridlen);
dY = reshape(dXYZ(2,:),gridlen,gridlen,gridlen);
dZ = reshape(dXYZ(3,:),gridlen,gridlen,gridlen);

%% Nonlinear transitions = nonunique vector field in R2
% dX = -ka*(X.*Y+X.*Z)        + kd*Z + fx*X;
% dY = -ka*(X.*Y+Y.*Z)        + kd*Z + fy*Y;
% dZ =  ka*(X.*Y+X.*Z+Y.*Z)   - kd*Z + fz*Z;

dN = dX+dY+dZ;
%Convert to fractional coordinates
dXf = (N.*dX - X.*dN)./(N.^2);
dYf = (N.*dY - Y.*dN)./(N.^2);
dZf = (N.*dZ - Z.*dN)./(N.^2);
figure
quiver3(Xf,Yf,Zf,dXf,dYf,dZf), view(135,0)
figure
mag = (dXf.^2+dYf.^2+dZf.^2).^(.5);
%Old -- colored scatter
%mag(isnan(mag))=0;
%scatter3(Xf(:),Yf(:),Zf(:),1000,mag(:),'filled'),view(135,0)
%mag(isnan(mag))=0;

%** Not all functions dXbar will yield unique vector fields on 2-simplex.
%   Some additionally depend on N (such as model above, hence multiple
%   arrows at each index)

%Simulate in 3-space, project back onto 2-simplex
T = [0 1 .5;0 0 1];

pts3 = [Xf(:)';Yf(:)';Zf(:)'];
pts2 = T*pts3;
tXf = pts2(1,:);
tYf = pts2(2,:);
pts3 = [dXf(:)';dYf(:)';dZf(:)'];
pts2 = T*pts3;
tdXf = pts2(1,:);
tdYf = pts2(2,:);

%New -- interpolate on normalized coordinates and plot heat (mesh)
mag(isnan(mag(:))|isinf(mag(:)))=0;
tXf(isnan(tXf)|isinf(tXf))=0;
tYf(isnan(tYf)|isinf(tYf))=0;
TSI = TriScatteredInterp(tXf',tYf',mag(:));
grid_fraction = linspace(0,1,1000);
[tXi,tYi] = meshgrid(grid_fraction, grid_fraction);
magi = TSI(tXi,tYi);
mesh(tXi,tYi,magi-max(magi(:))), view(0,90);
hold on;
quiver(tXf,tYf,tdXf,tdYf,'Color','k');

%Sample trajectory integration
[t,y] = ode23s(@(t,y) linear_payoff(t,y,A),[0 10],[1/3 1/3 1/3]);
yf = y./repmat(sum(y,2),1,3);
pts3 = yf'; 
tyf = T*pts3;
plot(tyf(1,:),tyf(2,:))

figure
plot(t,yf)
legend('x','y','z')

%add N as 3rd dimension
figure
quiver3(tXf,tYf,N(:)',tdXf,tdYf,dN(:)');

%% Demonstrate coordinate transformation
T = [0 1 .5;0 0 1];
[x,y] = meshgrid(0:.1:1,0:.1:1);
z=1-x-y;
x(z<0)=[];
y(z<0)=[];
z(z<0)=[];
pts3 = [x(:)';y(:)';z(:)'];
pts2 = [T*pts3;zeros(1,numel(x))];
figure
hold on
scatter3(pts3(1,:),pts3(2,:),pts3(3,:),'filled')
scatter3(pts2(1,:),pts2(2,:),pts2(3,:),'filled')
for i=1:size(pts3,2)
    line([pts3(1,i),pts2(1,i)],[pts3(2,i),pts2(2,i)],[pts3(3,i),pts2(3,i)],'Color',[.5 .5 .5]);
end
    
%%Symbolic
syms k fx fy fz A Xt X x y z xt(t) yt(t) zt(t);
%Linear payoffs, x+y->0+z, x+z->0+z, y+z->0+z
%dX = (A-phi)X, phi = Sum(dX)
    %x      %y    %z
%fy = fx;
A = [fx     -k     -k;    %x
     -k     fx     -k;    %y
     k      k      fz];   %z
X = [x;y;z];
Xt = [xt;yt;zt];
phi = ones(3,1)*transpose(X)*A;
[V,lambda] = eig(A);
%Y = dsolve(diff(Xt)==(A-phi)*Xt);
EQ = solve((A-phi)*X==0);
detA = det(A);
dX = (A-phi)*X;
J = jacobian((A-phi)*X,[x y z]);
EJ = eig(J);
EJn = subs(EJ,fx,1);
EJnfx = subs(EJn,[x y z],[subs(EQ.x(2),fx,1), subs(EQ.y(2),fx,1), subs(EQ.z(2),fx,1)]);
%Linear Stability - unclear
figure, hold on
ezplot(char(EJnfx(1,:)));
ezplot(char(EJnfx(2,:))); 
ezplot(char(EJnfx(3,:)));
axis([-2 2 0 2])
figure
ezcontourf(char(EJnfx(1,:),[0 5 0 5]));
figure
ezcontourf(char(EJnfx(2,:),[0 5 0 5]));
figure
ezcontourf(char(EJnfx(3,:),[0 5 0 5]));

%Parameter space for fixed points/no fixed points -- INCORRECT
%ezmesh('(fx^2 - 2*fx*fz - 2*fx*1 + fz^2 + 2*fz*1 - 7*1^2)>0',[0,10])
%a = k/fx, b=fz/fx
defaultPlot()
%ezcontourf('a - a*b + b',[0,2])
figure, hold on
ezplot('(fz + k)/(1 + 2*fz - k)')
ezplot('(fz + k)/(1 + 2*fz - k)-1')
ezplot('(1 - 3*k)/(1 + 2*fz - k)')
ezplot('(1 - 3*k)/(1 + 2*fz - k)-1')
%as lines
figure, hold on
line([-2 2],[-.5 1.5],'Color','r','LineWidth',2)
line([-2 0],[2,0],'Color','g','LineWidth',2)
line([-2 2],[1/3,1/3],'Color','b','LineWidth',2)
axis([-2 2 0 2])


%clarify
line([0 2],[1 1],'Color','k','LineWidth',2)
line([1 1],[0 2],'Color','k','LineWidth',2)
%title('k/f_x - k/f_x*f_z/f_x + f_z/f_x < 1')
xlabel('k/f_x')
ylabel('f_z/f_x')
set(gca(),'XTick',[0 1 2])
set(gca(),'YTick',[0 1 2])
colormap 'bone'
%Sample fixed points
subs([EQ.x EQ.y EQ.z]', [fx fz k], [7 3 1])
%No potential exists (has curl)
%potential(QV,X)

%manual map
a = 0:.1:2;
b = 0:.1:2;
fxpt = nan(size(a),size(b));
for i = 1:length(a)
    for j = 1:length(b)
        c = subs([EQ.x EQ.y EQ.z]', [fx fz k], [1 a(i) b(j)]);
        fxpt(i,j) = any(c(:)<0);
    end
end


QV = subs(dX,[fx fz k],[1 1 .5]);

dx = @(x,y,z) eval(vectorize(QV(1)));
dy = @(x,y,z) eval(vectorize(QV(2)));
dz = @(x,y,z) eval(vectorize(QV(3)));
grid_base = [0 logspace(0,2,10)];
[gx,gy,gz] = meshgrid(grid_base,grid_base,grid_base);
xpts = gx./(gx+gy+gz);
ypts = gy./(gx+gy+gz);
zpts = gz./(gx+gy+gz);
figure
quiver3(xpts,ypts,zpts,dx(gx,gy,gz),dy(gx,gy,gz),dz(gx,gy,gz));

%projection onto S3
Xf = xpts;
Yf = ypts;
Zf = zpts;
dXf = dx(gx,gy,gz);
dYf = dy(gx,gy,gz);
dZf = dz(gx,gy,gz);
pts3 = [Xf(:)';Yf(:)';Zf(:)'];
%--rewrite as function, but goto plotting step

AA = [          
    0, - fx - k, - fz - k
 - fx - k,          0, - fz - k
   k - fx,   k - fx,          0];



%% 4-state model (with 3 pure strategies)

